home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / POLYROOT.ZIP / PRQD.FOR < prev   
Text File  |  1985-11-29  |  10KB  |  388 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE PRQD
  5. C
  6. C        PURPOSE
  7. C           CALCULATE ALL REAL AND COMPLEX ROOTS OF A GIVEN POLYNOMIAL
  8. C           WITH REAL COEFFICIENTS.
  9. C
  10. C        USAGE
  11. C           CALL PRQD(C,IC,Q,E,POL,IR,IER)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           C     - COEFFICIENT VECTOR OF GIVEN POLYNOMIAL
  15. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  16. C                   THE GIVEN COEFFICIENT VECTOR GETS DIVIDED BY THE
  17. C                   LAST NONZERO TERM
  18. C           IC    - DIMENSION OF VECTOR C
  19. C           Q     - WORKING STORAGE OF DIMENSION IC
  20. C                   ON RETURN Q CONTAINS REAL PARTS OF ROOTS
  21. C           E     - WORKING STORAGE OF DIMENSION IC
  22. C                   ON RETURN E CONTAINS COMPLEX PARTS OF ROOTS
  23. C           POL   - WORKING STORAGE OF DIMENSION IC
  24. C                   ON RETURN POL CONTAINS THE COEFFICIENTS OF THE
  25. C                   POLYNOMIAL WITH CALCULATED ROOTS
  26. C                   THIS RESULTING COEFFICIENT VECTOR HAS DIMENSION IR+1
  27. C                   COEFFICIENTS ARE ORDERED FROM LOW TO HIGH
  28. C           IR    - NUMBER OF CALCULATED ROOTS
  29. C                   NORMALLY IR IS EQUAL TO DIMENSION IC MINUS ONE
  30. C           IER   - RESULTING ERROR PARAMETER. SEE REMARKS
  31. C
  32. C        REMARKS
  33. C           THE REAL PART OF THE ROOTS IS STORED IN Q(1) UP TO Q(IR)
  34. C           CORRESPONDING COMPLEX PARTS ARE STORED IN E(1) UP TO E(IR).
  35. C           IER = 0 MEANS NO ERRORS
  36. C           IER = 1 MEANS NO CONVERGENCE WITH FEASIBLE TOLERANCE
  37. C           IER = 2 MEANS POLYNOMIAL IS DEGENERATE (CONSTANT OR ZERO)
  38. C           IER = 3 MEANS SUBROUTINE WAS ABANDONED DUE TO ZERO DIVISOR
  39. C           IER = 4 MEANS THERE EXISTS NO S-FRACTION
  40. C           IER =-1 MEANS CALCULATED COEFFICIENT VECTOR REVEALS POOR
  41. C                   ACCURACY OF THE CALCULATED ROOTS.
  42. C                   THE CALCULATED COEFFICIENT VECTOR HAS LESS THAN
  43. C                   3 CORRECT DIGITS.
  44. C           THE FINAL COMPARISON BETWEEN GIVEN AND CALCULATED
  45. C           COEFFICIENT VECTOR IS PERFORMED ONLY IF ALL ROOTS HAVE BEEN
  46. C           CALCULATED.
  47. C           THE MAXIMAL RELATIVE ERROR OF THE COEFFICIENT VECTOR IS
  48. C           RECORDED IN Q(IR+1).
  49. C
  50. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  51. C           NONE
  52. C
  53. C        METHOD
  54. C           THE ROOTS OF THE POLYNOMIAL ARE CALCULATED BY MEANS OF
  55. C           THE QUOTIENT-DIFFERENCE ALGORITHM WITH DISPLACEMENT.
  56. C           REFERENCE
  57. C           H.RUTISHAUSER, DER QUOTIENTEN-DIFFERENZEN-ALGORITHMUS,
  58. C           BIRKHAEUSER, BASEL/STUTTGART, 1957.
  59. C
  60. C     ..................................................................
  61. C
  62.       SUBROUTINE PRQD(C,IC,Q,E,POL,IR,IER)
  63. C
  64. C      DIMENSIONED DUMMY VARIABLES
  65.       DIMENSION E(1),Q(1),C(1),POL(1)
  66. C
  67. C        NORMALIZATION OF GIVEN POLYNOMIAL
  68. C           TEST OF DIMENSION
  69. C        IR CONTAINS INDEX OF HIGHEST COEFFICIENT
  70.       IER=0
  71.       IR=IC
  72.       EPS=1.E-6
  73.       TOL=1.E-3
  74.       LIMIT=10*IC
  75.       KOUNT=0
  76.     1 IF(IR-1)79,79,2
  77. C
  78. C        DROP TRAILING ZERO COEFFICIENTS
  79.     2 IF(C(IR))4,3,4
  80.     3 IR=IR-1
  81.       GOTO 1
  82. C
  83. C           REARRANGEMENT OF GIVEN POLYNOMIAL
  84. C        EXTRACTION OF ZERO ROOTS
  85.     4 O=1./C(IR)
  86.       IEND=IR-1
  87.       ISTA=1
  88.       NSAV=IR+1
  89.       JBEG=1
  90. C
  91. C        Q(J)=1.
  92. C        Q(J+I)=C(IR-I)/C(IR)
  93. C        Q(IR)=C(J)/C(IR)
  94. C        WHERE J IS THE INDEX OF THE LOWEST NONZERO COEFFICIENT
  95.       DO 9 I=1,IR
  96.       J=NSAV-I
  97.       IF(C(I))7,5,7
  98.     5 GOTO(6,8),JBEG
  99.     6 NSAV=NSAV+1
  100.       Q(ISTA)=0.
  101.       E(ISTA)=0.
  102.       ISTA=ISTA+1
  103.       GOTO 9
  104.     7 JBEG=2
  105.     8 Q(J)=C(I)*O
  106.       C(I)=Q(J)
  107.     9 CONTINUE
  108. C
  109. C           INITIALIZATION
  110.       ESAV=0.
  111.       Q(ISTA)=0.
  112.    10 NSAV=IR
  113. C
  114. C        COMPUTATION OF DERIVATIVE
  115.       EXPT=IR-ISTA
  116.       E(ISTA)=EXPT
  117.       DO 11 I=ISTA,IEND
  118.       EXPT=EXPT-1.0
  119.       POL(I+1)=EPS*ABS(Q(I+1))+EPS
  120.    11 E(I+1)=Q(I+1)*EXPT
  121. C
  122. C        TEST OF REMAINING DIMENSION
  123.       IF(ISTA-IEND)12,20,60
  124.    12 JEND=IEND-1
  125. C
  126. C        COMPUTATION OF S-FRACTION
  127.       DO 19 I=ISTA,JEND
  128.       IF(I-ISTA)13,16,13
  129.    13 IF(ABS(E(I))-POL(I+1))14,14,16
  130. C
  131. C        THE GIVEN POLYNOMIAL HAS MULTIPLE ROOTS, THE COEFFICIENTS OF
  132. C        THE COMMON FACTOR ARE STORED FROM Q(NSAV) UP TO Q(IR)
  133.    14 NSAV=I
  134.       DO 15 K=I,JEND
  135.       IF(ABS(E(K))-POL(K+1))15,15,80
  136.    15 CONTINUE
  137.       GOTO 21
  138. C
  139. C           EUCLIDEAN ALGORITHM
  140.    16 DO 19 K=I,IEND
  141.       E(K+1)=E(K+1)/E(I)
  142.       Q(K+1)=E(K+1)-Q(K+1)
  143.       IF(K-I)18,17,18
  144. C
  145. C        TEST FOR SMALL DIVISOR
  146.    17 IF(ABS(Q(I+1))-POL(I+1))80,80,19
  147.    18 Q(K+1)=Q(K+1)/Q(I+1)
  148.       POL(K+1)=POL(K+1)/ABS(Q(I+1))
  149.       E(K)=Q(K+1)-E(K)
  150.    19 CONTINUE
  151.    20 Q(IR)=-Q(IR)
  152. C
  153. C           THE DISPLACEMENT EXPT IS SET TO 0 AUTOMATICALLY.
  154. C           E(ISTA)=0.,Q(ISTA+1),...,E(NSAV-1),Q(NSAV),E(NSAV)=0.,
  155. C           FORM A DIAGONAL OF THE QD-ARRAY.
  156. C        INITIALIZATION OF BOUNDARY VALUES
  157.    21 E(ISTA)=0.
  158.       NRAN=NSAV-1
  159.    22 E(NRAN+1)=0.
  160. C
  161. C           TEST FOR LINEAR OR CONSTANT FACTOR
  162. C        NRAN-ISTA IS DEGREE-1
  163.       IF(NRAN-ISTA)24,23,31
  164. C
  165. C        LINEAR FACTOR
  166.    23 Q(ISTA+1)=Q(ISTA+1)+EXPT
  167.       E(ISTA+1)=0.
  168. C
  169. C        TEST FOR UNFACTORED COMMON DIVISOR
  170.    24 E(ISTA)=ESAV
  171.       IF(IR-NSAV)60,60,25
  172. C
  173. C        INITIALIZE QD-ALGORITHM FOR COMMON DIVISOR
  174.    25 ISTA=NSAV
  175.       ESAV=E(ISTA)
  176.       GOTO 10
  177. C
  178. C        COMPUTATION OF ROOT PAIR
  179.    26 P=P+EXPT
  180. C
  181. C           TEST FOR REALITY
  182.       IF(O)27,28,28
  183. C
  184. C           COMPLEX ROOT PAIR
  185.    27 Q(NRAN)=P
  186.       Q(NRAN+1)=P
  187.       E(NRAN)=T
  188.       E(NRAN+1)=-T
  189.       GOTO 29
  190. C
  191. C           REAL ROOT PAIR
  192.    28 Q(NRAN)=P-T
  193.       Q(NRAN+1)=P+T
  194.       E(NRAN)=0.
  195. C
  196. C           REDUCTION OF DEGREE BY 2 (DEFLATION)
  197.    29 NRAN=NRAN-2
  198.       GOTO 22
  199. C
  200. C        COMPUTATION OF REAL ROOT
  201.    30 Q(NRAN+1)=EXPT+P
  202. C
  203. C           REDUCTION OF DEGREE BY 1 (DEFLATION)
  204.       NRAN=NRAN-1
  205.       GOTO 22
  206. C
  207. C        START QD-ITERATION
  208.    31 JBEG=ISTA+1
  209.       JEND=NRAN-1
  210.       TEPS=EPS
  211.       TDELT=1.E-2
  212.    32 KOUNT=KOUNT+1
  213.       P=Q(NRAN+1)
  214.       R=ABS(E(NRAN))
  215. C
  216. C           TEST FOR CONVERGENCE
  217.       IF(R-TEPS)30,30,33
  218.    33 S=ABS(E(JEND))
  219. C
  220. C        IS THERE A REAL ROOT NEXT
  221.       IF(S-R)38,38,34
  222. C
  223. C        IS DISPLACEMENT SMALL ENOUGH
  224.    34 IF(R-TDELT)36,35,35
  225.    35 P=0.
  226.    36 O=P
  227.       DO 37 J=JBEG,NRAN
  228.       Q(J)=Q(J)+E(J)-E(J-1)-O
  229. C
  230. C           TEST FOR SMALL DIVISOR
  231.       IF(ABS(Q(J))-POL(J))81,81,37
  232.    37 E(J)=Q(J+1)*E(J)/Q(J)
  233.       Q(NRAN+1)=-E(NRAN)+Q(NRAN+1)-O
  234.       GOTO 54
  235. C
  236. C        CALCULATE DISPLACEMENT FOR DOUBLE ROOTS
  237. C           QUADRATIC EQUATION FOR DOUBLE ROOTS
  238. C           X**2-(Q(NRAN)+Q(NRAN+1)+E(NRAN))*X+Q(NRAN)*Q(NRAN+1)=0
  239.    38 P=0.5*(Q(NRAN)+E(NRAN)+Q(NRAN+1))
  240.       O=P*P-Q(NRAN)*Q(NRAN+1)
  241.       T=SQRT(ABS(O))
  242. C
  243. C        TEST FOR CONVERGENCE
  244.       IF(S-TEPS)26,26,39
  245. C
  246. C        ARE THERE COMPLEX ROOTS
  247.    39 IF(O)43,40,40
  248.    40 IF(P)42,41,41
  249.    41 T=-T
  250.    42 P=P+T
  251.       R=S
  252.       GOTO 34
  253. C
  254. C        MODIFICATION FOR COMPLEX ROOTS
  255. C        IS DISPLACEMENT SMALL ENOUGH
  256.    43 IF(S-TDELT)44,35,35
  257. C
  258. C        INITIALIZATION
  259.    44 O=Q(JBEG)+E(JBEG)-P
  260. C
  261. C        TEST FOR SMALL DIVISOR
  262.       IF(ABS(O)-POL(JBEG))81,81,45
  263.    45 T=(T/O)**2
  264.       U=E(JBEG)*Q(JBEG+1)/(O*(1.+T))
  265.       V=O+U
  266.       KOUNT=KOUNT+2
  267. C
  268. C        THREEFOLD LOOP FOR COMPLEX DISPLACEMENT
  269.       DO 53 J=JBEG,NRAN
  270.       O=Q(J+1)+E(J+1)-U-P
  271. C
  272. C        TEST FOR SMALL DIVISOR
  273.       IF(ABS(V)-POL(J))46,46,49
  274.    46 IF(J-NRAN)81,47,81
  275.    47 EXPT=EXPT+P
  276.       IF(ABS(E(JEND))-TOL)48,48,81
  277.    48 P=0.5*(V+O-E(JEND))
  278.       O=P*P-(V-U)*(O-U*T-O*W*(1.+T)/Q(JEND))
  279.       T=SQRT(ABS(O))
  280.       GOTO 26
  281. C
  282. C           TEST FOR SMALL DIVISOR
  283.    49 IF(ABS(O)-POL(J+1))46,46,50
  284.    50 W=U*O/V
  285.       T=T*(V/O)**2
  286.       Q(J)=V+W-E(J-1)
  287.       U=0.
  288.       IF(J-NRAN)51,52,52
  289.    51 U=Q(J+2)*E(J+1)/(O*(1.+T))
  290.    52 V=O+U-W
  291. C
  292. C        TEST FOR SMALL DIVISOR
  293.       IF(ABS(Q(J))-POL(J))81,81,53
  294.    53 E(J)=W*V*(1.+T)/Q(J)
  295.       Q(NRAN+1)=V-E(NRAN)
  296.    54 EXPT=EXPT+P
  297.       TEPS=TEPS*1.1
  298.       TDELT=TDELT*1.1
  299.       IF(KOUNT-LIMIT)32,55,55
  300. C
  301. C        NO CONVERGENCE WITH FEASIBLE TOLERANCE
  302. C           ERROR RETURN IN CASE OF UNSATISFACTORY CONVERGENCE
  303.    55 IER=1
  304. C
  305. C        REARRANGE CALCULATED ROOTS
  306.    56 IEND=NSAV-NRAN-1
  307.       E(ISTA)=ESAV
  308.       IF(IEND)59,59,57
  309.    57 DO 58 I=1,IEND
  310.       J=ISTA+I
  311.       K=NRAN+1+I
  312.       E(J)=E(K)
  313.    58 Q(J)=Q(K)
  314.    59 IR=ISTA+IEND
  315. C
  316. C        NORMAL RETURN
  317.    60 IR=IR-1
  318.       IF(IR)78,78,61
  319. C
  320. C        REARRANGE CALCULATED ROOTS
  321.    61 DO 62 I=1,IR
  322.       Q(I)=Q(I+1)
  323.    62 E(I)=E(I+1)
  324. C
  325. C        CALCULATE COEFFICIENT VECTOR FROM ROOTS
  326.       POL(IR+1)=1.
  327.       IEND=IR-1
  328.       JBEG=1
  329.       DO 69 J=1,IR
  330.       ISTA=IR+1-J
  331.       O=0.
  332.       P=Q(ISTA)
  333.       T=E(ISTA)
  334.       IF(T)65,63,65
  335. C
  336. C        MULTIPLY WITH LINEAR FACTOR
  337.    63 DO 64 I=ISTA,IR
  338.       POL(I)=O-P*POL(I+1)
  339.    64 O=POL(I+1)
  340.       GOTO 69
  341.    65 GOTO(66,67),JBEG
  342.    66 JBEG=2
  343.       POL(ISTA)=0.
  344.       GOTO 69
  345. C
  346. C        MULTIPLY WITH QUADRATIC FACTOR
  347.    67 JBEG=1
  348.       U=P*P+T*T
  349.       P=P+P
  350.       DO 68 I=ISTA,IEND
  351.       POL(I)=O-P*POL(I+1)+U*POL(I+2)
  352.    68 O=POL(I+1)
  353.       POL(IR)=O-P
  354.    69 CONTINUE
  355.       IF(IER)78,70,78
  356. C
  357. C        COMPARISON OF COEFFICIENT VECTORS, IE. TEST OF ACCURACY
  358.    70 P=0.
  359.       DO 75 I=1,IR
  360.       IF(C(I))72,71,72
  361.    71 O=ABS(POL(I))
  362.       GOTO 73
  363.    72 O=ABS((POL(I)-C(I))/C(I))
  364.    73 IF(P-O)74,75,75
  365.    74 P=O
  366.    75 CONTINUE
  367.       IF(P-TOL)77,76,76
  368.    76 IER=-1
  369.    77 Q(IR+1)=P
  370.       E(IR+1)=0.
  371.    78 RETURN
  372. C
  373. C        ERROR RETURNS
  374. C           ERROR RETURN FOR POLYNOMIALS OF DEGREE LESS THAN 1
  375.    79 IER=2
  376.       IR=0
  377.       RETURN
  378. C
  379. C           ERROR RETURN IF THERE EXISTS NO S-FRACTION
  380.    80 IER=4
  381.       IR=ISTA
  382.       GOTO 60
  383. C
  384. C           ERROR RETURN IN CASE OF INSTABLE QD-ALGORITHM
  385.    81 IER=3
  386.       GOTO 56
  387.       END
  388.